library(sf)
library(stringr)
library(dplyr)
library(sp)
library(raster)
library(rgeos)
library(rgdal)
library(rasterVis)
library(ggplot2)
library(terra)
library(sf)
library(fs)
library(here)
library(dplyr)
# TEST IMAGES
image1B5 <- raster(here("landsat_images","LC08_L2SP_027041_20180829_20200831_02_T1","LC08_L2SP_027041_20180829_20200831_02_T1_SR_B5.TIF"))
image1B4 <- raster(here("landsat_images","LC08_L2SP_027041_20180829_20200831_02_T1","LC08_L2SP_027041_20180829_20200831_02_T1_SR_B4.TIF"))
image1B2 <- raster(here("landsat_images","LC08_L2SP_027041_20180829_20200831_02_T1","LC08_L2SP_027041_20180829_20200831_02_T1_SR_B2.TIF"))
image1B6 <- raster(here("landsat_images","LC08_L2SP_027041_20180829_20200831_02_T1","LC08_L2SP_027041_20180829_20200831_02_T1_SR_B6.TIF"))
image2B5 <- raster(here("landsat_images","LC08_L2SP_027041_20170826_20200903_02_T1","LC08_L2SP_027041_20170826_20200903_02_T1_SR_B5.TIF"))
image2B4 <- raster(here("landsat_images","LC08_L2SP_027041_20170826_20200903_02_T1","LC08_L2SP_027041_20170826_20200903_02_T1_SR_B4.TIF"))
image2B2 <- raster(here("landsat_images","LC08_L2SP_027041_20170826_20200903_02_T1","LC08_L2SP_027041_20170826_20200903_02_T1_SR_B2.TIF"))
image2B6 <- raster(here("landsat_images","LC08_L2SP_027041_20170826_20200903_02_T1","LC08_L2SP_027041_20170826_20200903_02_T1_SR_B6.TIF"))
Pre1B5 <- raster(here("landsat_images","LC08_L2SP_027041_20170911_20200903_02_T1","LC08_L2SP_027041_20170911_20200903_02_T1_SR_B5.TIF"))
Pre1B4 <- raster(here("landsat_images","LC08_L2SP_027041_20170911_20200903_02_T1","LC08_L2SP_027041_20170911_20200903_02_T1_SR_B4.TIF"))
Pre1B2 <- raster(here("landsat_images","LC08_L2SP_027041_20170911_20200903_02_T1","LC08_L2SP_027041_20170911_20200903_02_T1_SR_B2.TIF"))
Pre1B6 <- raster(here("landsat_images","LC08_L2SP_027041_20170911_20200903_02_T1","LC08_L2SP_027041_20170911_20200903_02_T1_SR_B6.TIF"))
#RGB_image <- stack(image1B5,image1B4,image1B2)
#RGB_image <- clamp(RGB_image,0,65535)
#RGB_image <-floor((RGB_image/65535)*255)
#plotRGB(RGB_image,axes=TRUE, main = "False Color: 5,4,2" )
RGB_image2 <- stack(image2B5,image2B4,image2B2)
RGB_image2 <- clamp(RGB_image2,0,65535)
RGB_image2 <-floor((RGB_image2/65535)*255)
plotRGB(RGB_image2,axes=TRUE, main = "False Color: 5,4,2" )

RGB_image <- stack(image1B5,image1B4,image1B2)
RGB_image <- clamp(RGB_image,0,65535)
RGB_image <-floor((RGB_image/65535)*255)
plotRGB(RGB_image,axes=TRUE, main = "False Color Post Hurricane 5,4,2" )

Pre1RGB <- stack(Pre1B5,Pre1B4,Pre1B2)
Pre1RGB <- clamp(Pre1RGB,0,65535)
Pre1RGB <- floor((Pre1RGB/65535)*255)
plotRGB(Pre1RGB,axes=TRUE, main="Fasle Color Pre Hurricane 5,4,2")

#NDVI and NDII analysis
NDVIfun <- function(NIR, red) {
NDVI <- (NIR - red) / (NIR + red)
return(NDVI)
}
# Uses SWIR band instead of NIR to look at water absoprtion levels and draw a contrast between them to distinguish rich vegetation
NDIIfun <- function(SWIR,red) {
NDII <- (SWIR - red) / (SWIR + red)
return(NDII)
}
#NDII
ndvi <- NDVIfun(Pre1B5,Pre1B4)
ndvi %>%
plot(.,col=rev(terrain.colors(10)),main ="Landsat NDVI Pre Hurricane Harvey")

NdviPost <- NDVIfun(image1B5,image1B4)
NdviPost %>%
plot(.,col=rev(terrain.colors(10)),main ="Landsat NDVI Post Hurricane Harvey")

ndvi %>%
hist(.,breaks = 40, main="NDVI Pre Histogram",xlim =c(-.3,.8))
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

NdviPost%>%
hist(.,breaks=40,main="NDVI Post Histogram",xlim=c(-.3,.8))
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

#NDII
ndiipre <- NDIIfun(Pre1B6,Pre1B4)
ndiipre %>%
plot(.,col=rev(terrain.colors(10)),main ="Landsat NDII Pre Hurricane Harvey")

ndiipost <- NDIIfun(image1B6,image1B4)
ndiipost %>%
plot(.,col=rev(terrain.colors(10)),main ="Landsat NDII Post Hurricane Harvey")

ndiipre %>%
hist(.,breaks=40,main="NDII Pre Histogram",xlim=c(-.3,.8))
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

ndiipost %>%
hist(.,breaks=40,main="NDII Post Histogram",xlim=c(-.3,.8))

PreMarB5 <- raster(here("landsat_images","LC08_L2SP_005047_20171003_20200903_02_T1","LC08_L2SP_005047_20171003_20200903_02_T1_SR_B5.TIF"))
PreMarB4 <- raster(here("landsat_images","LC08_L2SP_005047_20171003_20200903_02_T1","LC08_L2SP_005047_20171003_20200903_02_T1_SR_B4.TIF"))
PreMarB2 <- raster(here("landsat_images","LC08_L2SP_005047_20171003_20200903_02_T1","LC08_L2SP_005047_20171003_20200903_02_T1_SR_B2.TIF"))
PreMarB6<-raster(here("landsat_images","LC08_L2SP_005047_20171003_20200903_02_T1","LC08_L2SP_005047_20171003_20200903_02_T1_SR_B6.TIF"))
PostMarB5 <- raster(here("landsat_images","LC08_L2SP_005047_20180904_20200831_02_T1","LC08_L2SP_005047_20180904_20200831_02_T1_SR_B5.TIF"))
PostMarB4 <- raster(here("landsat_images","LC08_L2SP_005047_20180904_20200831_02_T1","LC08_L2SP_005047_20180904_20200831_02_T1_SR_B4.TIF"))
PostMarB2 <- raster(here("landsat_images","LC08_L2SP_005047_20180904_20200831_02_T1","LC08_L2SP_005047_20180904_20200831_02_T1_SR_B2.TIF"))
PostMarB6 <- raster(here("landsat_images","LC08_L2SP_005047_20180904_20200831_02_T1","LC08_L2SP_005047_20180904_20200831_02_T1_SR_B6.TIF"))
Mar1 <- stack(PreMarB5, PreMarB4,PreMarB2)
Mar1 <- clamp(Mar1,0,65535)
Mar1 <-floor((Mar1/65535)*255)
plotRGB(Mar1,axes=TRUE, main = "False Color Pre Maria Hurricane 5,4,2" )

Mar2 <- stack(PostMarB5,PostMarB4,PostMarB2)
Mar2 <- clamp(Mar2,0,65535)
Mar2 <- floor((Mar2/65535)*255)
plotRGB(Mar2,axes=TRUE,main="False Color Post Maria Hurricane 5,4,2")

ndvipreMar <- NDVIfun(PreMarB5,PreMarB4)
ndvipreMar %>%
plot(.,col=rev(terrain.colors(10)),main ="Landsat NDVI Pre Hurricane Maria")

ndvipostMar <- NDVIfun(PostMarB5,PostMarB4)
ndvipostMar %>%
plot(.,col=rev(terrain.colors(10)),main="Landsat NDVI Post Hurricane Maria")

ndvipostMar %>%
hist(.,breaks=40,main="Histrogram NDVI Post Maria")

ndvipreMar%>%
hist(.,breaks=40,main="Histogram NDVI Pre Maria")

ndiipremar <- NDIIfun(PreMarB6,PreMarB4)
ndiipremar %>%
plot(.,col=rev(terrain.colors(10)),main="Landsat NDII Pre Hurricane Maria")

ndiipostmar <- NDIIfun(PostMarB6,PostMarB4)
ndiipostmar%>%
plot(.,col=rev(terrain.colors(10)),main="Landsat NDII Post Hurricane Maria")

ndiipostmar %>%
hist(.,breaks=40,main="Histogram NDII Pre Maria")
## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 0%
## of the raster cells were used. 100000 values used.

ndiipremar%>%
hist(.,breaks=40,main="Histrogram NDII Post Maria")

PreIan5 <- raster(here("landsat_images","LC08_L2SP_016041_20211011_20211019_02_T1","LC08_L2SP_016041_20211011_20211019_02_T1_SR_B5.TIF"))
PreIan4 <- raster(here("landsat_images","LC08_L2SP_016041_20211011_20211019_02_T1","LC08_L2SP_016041_20211011_20211019_02_T1_SR_B4.TIF"))
PreIan2 <- raster(here("landsat_images","LC08_L2SP_016041_20211011_20211019_02_T1","LC08_L2SP_016041_20211011_20211019_02_T1_SR_B2.TIF"))
PreIan6 <- raster(here("landsat_images","LC08_L2SP_016041_20211011_20211019_02_T1","LC08_L2SP_016041_20211011_20211019_02_T1_SR_B6.TIF"))
PostIan5 <- raster(here("landsat_images","LC09_L2SP_016041_20221107_20230323_02_T1","LC09_L2SP_016041_20221107_20230323_02_T1_SR_B5.TIF"))
PostIan4 <- raster(here("landsat_images","LC09_L2SP_016041_20221107_20230323_02_T1","LC09_L2SP_016041_20221107_20230323_02_T1_SR_B4.TIF"))
PostIan2 <- raster(here("landsat_images","LC09_L2SP_016041_20221107_20230323_02_T1","LC09_L2SP_016041_20221107_20230323_02_T1_SR_B2.TIF"))
PostIan6 <- raster(here("landsat_images","LC09_L2SP_016041_20221107_20230323_02_T1","LC09_L2SP_016041_20221107_20230323_02_T1_SR_B6.TIF"))
Ian1 <- stack(PreIan5, PreIan4,PreIan2)
Ian1 <- clamp(Ian1,0,65535)
Ian1 <-floor((Ian1/65535)*255)
plotRGB(Ian1,axes=TRUE, main = "False Color Pre Ian Hurricane 5,4,2" )

Ian2 <- stack(PostIan5,PostIan4,PostIan2)
Ian2 <- clamp(Ian2,0,65535)
Ian2 <- floor((Ian2/65535)*255)
plotRGB(Ian2,axes=TRUE,main="False Color Post Ian Hurricane 5,4,2")

ndvipreIAN <- NDVIfun(PreIan5,PreIan4)
ndvipreIAN%>%
plot(.,col=rev(terrain.colors(10)),main="NDVI Pre Hurricane Ian")

ndvipreIAN%>%
hist(.,breaks=40,main="Histogram NDVI Pre Ian")

ndvipostIAN <- NDVIfun(PostIan5,PostIan4)
ndvipostIAN%>%
plot(.,col=rev(terrain.colors(10)),main="NDIV Post Hurricane Ian")

ndvipostIAN%>%
hist(.,breaks=40,main="Histogram NDVI Post IAN")

ndiipreIAN <- NDIIfun(PreIan6,PreIan4)
ndiipreIAN %>%
plot(.,col=rev(terrain.colors(10)),main="NDII Pre Hurricane Ian")

ndiipreIAN%>%
hist(.,breaks=40,main="Histogram NDII PRE IAN")

ndiipostIan <- NDIIfun(PostIan6,PostIan4)
ndiipostIan%>%
plot(.,col=rev(terrain.colors(10)),main="NDII Post Hurricane Ian")

ndiipostIan%>%
hist(.,breaks=40,main="Histogram NDII POST IAN")

# magian <- ndiipostIan - ndiipreIAN
#magian %>%
# plot(.,col=rev(terrain.colors(20)))
magian2 <- ndvipostIAN - ndvipreIAN
## Warning in ndvipostIAN - ndvipreIAN: Raster objects have different extents.
## Result for their intersection is returned
#magian2%>%
# plot(.,col=rev(terrain.colors(20)))
layout_matrix <- matrix(c(1, 2, 3, 4,5,6), nrow = 3, ncol = 2, byrow = TRUE)
layout(layout_matrix)
magmar <- ndvipreMar - ndvipostMar
## Warning in ndvipreMar - ndvipostMar: Raster objects have different extents.
## Result for their intersection is returned
maghar <- NdviPost - ndvi
magianii <- ndiipostIan- ndiipreIAN
## Warning in ndiipostIan - ndiipreIAN: Raster objects have different extents.
## Result for their intersection is returned
magmarii <- ndiipostmar - ndiipremar
## Warning in ndiipostmar - ndiipremar: Raster objects have different extents.
## Result for their intersection is returned
magharii <- ndiipre - ndiipost
hist(magian2, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDVI Difference", ylab = "Number of Pixels", main = "Difference in NDVI")
hist(magmar, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDVI Difference", ylab = "Number of Pixels", main = "Difference in NDVI")
hist(maghar, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDVI Difference", ylab = "Number of Pixels", main = "Difference in NDVI")
hist(magianii, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDII Difference", ylab = "Number of Pixels", main = "Difference in NDII")
hist(magmarii, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDII Difference", ylab = "Number of Pixels", main = "Difference in NDII")
hist(magharii, breaks = 40,xlim=c(-0.4,0.6), xlab = "NDII Difference", ylab = "Number of Pixels", main = "Difference in NDII")
